Field dependence of the vortex structure in chiral p-wave superconductors 
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To investigate the different vortex structure between two chiral pairing p x ± ip y , we calculate the 
pair potential, the internal field, the local density of states, and free energy in the vortex lattice 
state based on the quasiclassical Eilenberger theory, and analyze the magnetic field dependence. 
The induced opposite chiral component of the pair potential plays an important role in the vortex 
structure. It also produces \/i?-behavior of the zero-energy density of states at higher field. These 
results are helpful when we understand the vortex states in Sr2Ru04. 
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I. INTRODUCTION 

For the superconducting state in quasi-two- 
dimensional metal Si^RuO.^, the pairing symmetry 
is suggested to be the chiral p-wave pairing with the 
basic-form A± ~ p x ± ip y and inplane equal-spin pair- 
ing.B@ For the experimental evidence, the spin triplet 
pairing is supported by the 17 0-NMR measurement, 
which reported that there is no reduction of the Knight 
shift in the superconducting state.ta The pairing state 
with broken time reversal symmetry is claimed by the 
/iSR measurement, which reported that spontaneous 
magnetic moment is induced in the superconducting 
state.Q 

Since the A + state and the A_ state are degenerate 
at zero field, we expect the domain formation of the two 
states A + and A_. We call them as the p + -wave domain 
and the p_-wave domain, respectively. This degeneracy 
is lifted under external magnetic field perpendicular to 
the basal plane, since A± is the broken time reversal 
symmetry state with an orbital angular momentum along 
the z-axis. Then, the vortex in the mixed state shows 
the different structure for the p + -wave and the p_-wave 
domains. We consider the case when H || z for e > (2e 
is the charge of the Cooper pair) or equivalently when 
H || — z for e < 0. In these cases, there appear winding 
+1 vortices. It is parallel (antiparallel) to the internal 
winding of the Cooper pair in the p + -wave (p_-wave) 
domain. The information of the vortex structure for the 
p ± -wave domain is important to analyze the chirality p± 
of each domain. 

The differences of the vortex structure for the p+-wave 
and the p_-wave domains have been studictLby the two 
component Ginabprg-Laudau (GL) theoryjju the quasi- 
classicaLthfiprylla and the Bogoliubov-de Gennes (BdG) 
theory.ffllilJ'El In the chiral p-wave superconductors, it is 
important to consider the opposite chiral component A T 
which is induced around the vortex of the dominant A± 
component in the p±-w&ve domains. This induced com- 
ponent shows different spatial structure for the p+-wave 
and the p_-wave domains. This is the origin of the dif- 
ferent vortex structure. There is a free energy difference 
between the p+-wave and the p_-wave domain cases in 
the vortex state, which leads to the different upper criti- 



cal field H C 2 for these chiral states. The estimation of H C 2 
was reported by Scharnberg and Kle mm .in the isotropic 
three dimensional Fermi surface case.tZlliia Following their 
results, H C 2 of the p + -wave domain, which is noted as 
ABM state, is near that of the isotropic s-wave pairing 
case, because the opposite A_-component is not induced 
along Hc2-line. On the other hand, the vortex state in 
the p_-wave domain, which is noted as the generalized 
ABM state or the SK state, gives higher H c2 because 
there is large induced A + -componcnt. That is, we can 
say that the superconductivity survives until high field 
in the p_-wave domain case by the enhancement due to 
the induced opposite chiral component. Since the differ- 
ence of H C 2 is large (about twice) for these domain cases, 
it is important to know the H (external magnetic field) 
dependence of the vortex structure for A + and A_ con- 
tinuously for all H regions, in order to study the chiral- 
dependent properties. 

In this paper, we investigate the field dependence of 
the vortex structure for the p±-wave domain cases, based 
on the quasiclassical Eilenberger theoryO Our calcula- 
tion method for the vortex lattice state was established 

for the study of ^f^rfr-w^i^Wf P arrm g case m high- 
T c superconductors .L300£§I1jIH3 We can calculate both 
the pair potential and the vector potential sclfconsis- 
tently. In this paper, we consider the case-of small GL 
parameter n = 2.7 appropriate to S^RuO^lj In the qua- 
siclassical theory, we can also consider the quasiparti- 
cle states around the vortex, such as the local density 
of states (LDOS). Among them, iJ-dependence of the 
spatially averaged zero energy density of states (DOS) 
iV(0) is important. It is accessible by the specific heat 
measurement. For example, in the d x 2_ y 2-w&ve pairing 
such as in high T c superconductors, there is a relation 
N(0) ~ spH. due._ta_the line node of the superconducting 
gap for H || C.E3E30 It is because low energy quasipar- 
ticles propagating to the node direction can extend far 
from vortex. In the relation N(0) ~ H a in the s-wave 
pairing case, a can be smaller than 1. It isj£la±ed_ta 
the field dependence of the vortex core radius .E§EjE3B3 E£I 
These behaviors could be confirmed by the calculation 
based on the quasiclassical theory.113 In Sr2Ru04, the spe- 
cific heat measurement suggests the relation N(0) ~ VH 
at higher field, giving the discussion on the possibility of 
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the line node@ Then, it is important to examine the 
origin of the \/fi-behavior. 

There are detailed discussions on the pairing func- 
tion of the chiral p-wave for Sr2Ru04, iiael 
anisotropy and the orbital dependence^' 
But, here, we consider the fundamental case of the sim- 
ple isotropic gap function A± = (p x ± ip y )/pY = e ±l9 on 
the two-dimensional isotropic Fermi surface (p x ,p y ) = 
Pf (cos 9, sin 8) in order to focus on the chirality depen- 
dence, without including the anisotropy effect of the su- 
perconducting gap or Fermi surfaces. 

After describing our formulation of the quasiclassical 
theory in Sec. II, we evaluate the free energy in Sec. [II. 



given as 



The pair potential and internal field structure of the vor- 
tex lattice state are studied in Sec. |y|. The low energy 
quasiparticle states are examined in Sec. by consider- 
ing the LDOS and the field dependence of the DOS. The 
last section is devoted to summary and discussions. 



II. QUASICLASSICAL EILENBERGER 
THEORY 



Our calculation is performed by extending the quasi- 
classical method fmJiix^yQrtey lattice state in the d x 2_ y 2- 
wave pairing casai3E30li3'E3a to the chiral p-wave pair- 
ing case. For the details of the calculation method, 
also see Refs. ^,^.^0[ We consider the case of the 
clean limit and cylindrical Fermi surface. First, to ob- 
tain the pair potential A(6*,r) and vector potential A(r) 
selfconsistently, we solve the Eilenberger equation in the 
Matsubara frequency uj n — (2n + 1)ttT for the quasi- 
classical Green's functions g(iuj n ,9, r), f(iui n ,9,r) and 
f'(ioj n , 9, r), where r is the center of mass coordinate of 
a Cooper pair. The direction of the relative momentum 
of the Cooper pair, k = k/|k|, is denoted by an angle 9 
measured-from the x axis. The Eilenberger equation is 
given byo 



^ F -(? + ^A(r))}/(io, n ,0,r) 



i. 2 v i <po 
= A(0,r)g(]w n ,6,r), 



(1) 



I 2 \ i <p 
= A*(9,r)g(iiu n ,9,v), (2) 
g(iw n , 9, r) = [1 - /(iw„, 9, r)/t(iu>„, 9, r)] 1 ^ ( 3 ) 

where Keg(iu n , 9, r) > 0, the Fermi velocity vp = vpk., 
the flux quantum </>o and e = — |e|. In the symmetric 
gauge, A(r) = x r + a(r), where H = (0,0, If) is 
a uniform field and a(r) is related to the internal field 
h(r) = (0,0, h(r)) as h(r) = V x a(r). For the coupling 
to a magnetic field, we neglect the paramagnetic coupling 
to the spin, and consider the effect of the orbital coupling 
in the vector potential terms. 

The self-consistent conditions for A(9, r) and a(r) are 



A(0,r) = N 2ifT V f 



2^ 



V (9', 9) j ; (iw n ,0',r), (4) 



V x V x a(r) 



d8k 
2ttT 



g(iu> n ,9,r) 



(5) 

with the pairing interaction V(6',6) and n = 
(7C(3)/72) 1 / 2 (A /T c )k BC s with Rieman's zeta function 
C(3). No is the density of states at the Fermi surface, Ao 
is the uniform gap at T — 0, and kbcs is the GL parame- 
ter in the BCS theory. We set the energy cutoff u c = 20T C 
and kbcs = 2.7. In the following, energies and lengths 
are measured in units of Ao and £o = vf/Aq = 7t£bcs 
(£bcs is the BCS coherence length), respectively. The 
magnetic fields are measured in units of </>o/£o- 

By solving Eqs. (|l|)-(§) in the so-called explosion 
method, we estimate the quasiclassical Green's functions 
at 41 x 41 discretized points in a unit <pUr-Qf the vor- 
tex lattice. Using the symmetry relationE3E3 described 
in Appendix, we can reduce the range of r and 9 in the 
calculation solving Eqs. We obtain new A(6, r) 

and a(r) from Eqs. (Q) and (p|), and use them at the 
next step calculation of Eqs. (|l])-(||). This iteration pro- 
cedure is repeated until sufficiently selfconsistent solution 
is obtained. When we consider the lattice transformation 
R = mri + 77X2 (to, n : integers) with the unit vectors 
ri = (a x ,0), ri = (C,a x ,a y ) of the vortex lattice and 
Ha x a y — 4>o, there is a relation 

A(0,r + R) = A(0,r)e ix(r ' R) , a(r + R)=a(r), (6) 
where 

x(r,R) = --^(HxR)-(r + 2r )-^TOn (7) 

in the symmetric gauge. Then, we can know A(9, r) and 
a(r) in the other region out of the calculated unit cell re- 
gion. There is a vortex center at ro — 5 (ri+r2). When we 
consider the case when a vortex center locates at r = 0, 
we set ro — i(ri + r2). The spatial variation of the in- 
ternal field and the current J(r) = (c/47r)V x h(r) is 
calculated from a(r). The current is scaled by c</>o/47r£ 
in figures. To study the field dependence, our calcula- 
tions are performed for various fields at fixed tempera- 
ture T/T c = 0.5. 

In the chiral p-wave pairing, we can set 



A(M)=A+(r)^+(0)- 

v{9',e)=v[4>* + {e')cj> + {e) 



A_(r) 



-(*)> 



(8) 
(9) 



with the pairing functions 4>±{9) — e ±w . For the p + - 
wave domain case, we start our calculations from the 
initial state A + (r) = \&o( r ) an d A_(r) = with the vor- 
tex lattice solution ^(r) in the lowest Landau level. In 
this case, A_(r) gives the induced component around 
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the vortex after we obtain selfconsistent results. For 
the p_-wave domain case, we start from the initial state 
A_(r) = * (r) and A+(r) = 0. 

When we consider the one component case for the 
pair potential by neglecting the induced component, 
A((9,r) = A+(r)</> + (0) and A(0,r) = A_(r)0_(0) in the 
p + -wave and the p_-wave domain cases, respectively. In 
this case, by setting <f>±{6) — \4>±(&)\e ±ie , /(iw n ,0,r) = 
f(icj n ,8,r)e ±ie and /t(iw n ,0,r) = pQu n , 0, r)e^, we 
can remove the phase factor e ±I0 out of the Eilen- 
bcrgcr equations. Then, the Eilenberger equations for 
/(iw ra , 0, r) and /t(iw ra , 0, r) are solved under the pair po- 
tential A(0,r) = A±(r)|0±(0)| = A±(r), which is the 
pair potential for the isotropic s-wave pairing. Then, 
unless we consider the induced component, the vortex 
structure is the same as that of the s-wave pairing case, 
and there are no differences for the p + -wave and the p_- 
wave domain cases. Therefore, two component pair po- 
tential is intrinsic and essential for the vortex structure 
in the chiral p-wave superconductors. We also calculate 
the vortex structure in the isotropic s-wave pairing case 
for reference. 

Tlxe free energy F is calculated asoE3L3 



2 (h(r) 2 )r 

(<t> /e ) 2 

"AT 

^° a>»>0 



|A + (r)| 2 + |A_(r)| 5 



2 71 



d0 
2^ 



N VA 2 
(JK,0,r)) r 



(10) 



with 



l(u n ,0,r) =A*( 

1 



7^ 

A*(0,r)/ + A(0,r)/t 



,r)/ + A(0,r)/t 
V 2tt 



2 VF 



(ii) 

(12) 



where <?, / and f' mean g(ioj n ,0,r), /(iw n ,0, r) and 
P(iui n , 0, r), respectively. For the spatial average, 
(• • •), = f ,, dr(- • -)/S, where S is the area of a unit 

\ > r J unit cell X, JJ. 1 . , 

cell. We use Eqs. (|l|)-(g|) to obtain Eq. ([12]). 
The LDOS for energy E is calculated as 



N(E,r) =N Q 



d0 
2^ 



Re g(iuj n E + ir),6,r), (13) 



To obtain g(iui n — > E + iry, 0, r), we solve Eqs. (|l|)-(^) for 
X] — IE instead of uj n using the self-consistently obtained 
A(0,r) and a(r). We typically use tj = 0.01. The DOS 
is given by the spatial average of the LDOS as 



N(E) = (N(E,T)) t 



(14) 



-o.i - 



-0.2 



H 



FIG. 1: Free energy as a function of H. We plot F/(N A%) 
for the p+-wave domain case (o) and the p_-wave domain case 
(•) • 



III. FREE ENERGY 

For the estimation of the free energy difference between 
the p + -wave domain and the p_-wave domain cases, we 
show the field dependence of the free energy F in Fig. 
0. The p_-wave domain case has lower free energy than 
that of the p + -wave domain case. Then, the p_-wave do- 
main is stable, and the p+-wave domain may exist as a 
metastable state. We expect that the transition of the 
p+-wave domain to the p_-wave domain is stimulated 
with increasing field, since the free energy difference in- 
creases. The upper critical field H c2 ~ 0.81 in the p + - 
wave domain case, and H c2 ~ 1.7 in the p_-wave do- 
main case. Then, the p+-wave domain does not exist at 
H > 0.81. 

The estimation of the stable vortex lattice configu- 
ration is also important, since the vortex lattice may 
be deformed from the triangular, lattice due to the in- 
duced opposite chiral component £3li3 Since our calcula- 
tion method needs long computational time, we can not 
check all possible vortex lattice configurations. Then, 
we compare the free energy for the triangular and the 
square lattice configurations, and discuss the stable vor- 
tex lattice configuration. The free energy difference is 
less than 10~ 3 of F. We use finer 121 x 121 mesh within 
a unit cell to carefully estimate the difference. For higher 
field H > 0.35, the square lattice configuration has lower 
free energy than the triangular one. It suggests that the 
square lattice is stable at higher field in the p_-wave do- 
main case. This is qualitatively consistent to the analy- 
sis by the GL theory, which-suggests the square lattice 
configuration at higher field.Oc3 And the square lattice 
is observed by neutron scattering experiment .c3 On the 
other hand, we obtain the result that the free energy of 
the triangular lattice is lower for all H range in the p+- 
wave domain case. It is because the induced opposite 
chiral component is small, as discussed later, and the 
vortex structure is similar to that of the isotropic s-wave 
pairing case. If the p_-wave domain and the metastable 
p+-wave domain coexist, we may observe the different 
vortex lattice configurations for the domains. 

In the following subsections, we investigate the ori- 
gin of the difference of the vortex structure between the 
p_-wave domain and the p+-wave domains. The vortex 
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structure is examined both in the square and the tri- 
angular vortex lattice configurations. Since our purpose 
is to find the chirality effect on the vortex structure by 
comparing the both chirality cases, we mainly report the 
results in the same situation of the square lattice con- 
figuration. The square lattice is expected in the stable 
p_-wave domain case. After that, we briefly comment on 
the triangular lattice configuration case. 



IV. VORTEX STRUCTURE 
A. Vortex structure in the p -wave domain 

First, we consider the vortex structure in the p_-wave 
domain case. It is shown in Fig. || (a) and (b) within 
a unit cell of the square vortex lattice, i.e., square area 
in Fig. ^ (a). The profiles are presented along the path 
V-S-C-V shown in Fig. || (a). Line VS is along the near- 
est neighbor (NN) vortex direction, and line SC is along 
the boundary of the unit cell. Dashed line VC is along 
the next nearest neighbor (NNN) vortex direction. The 
dominant component |A_(r)| shows a conventional vor- 
tex structure. At low field, as shown in Fig. ||(a), | A_ (r) | 
is recovered to Ao outside of the vortex core. The shape 
of the vortex core is square-like. At higher field H = 0.8 
(~ 0.5iJ C 2) shown in Fig. |^(b), |A_(r)| is not recovered 
to Ao even in the boundary of the unit cell, since the 
inter-vortex distance is small. Along the NN vortex di- 
rection, |A_(r)| is slightly suppressed at the boundary 
region, compared to the NNN vortex direction. 

The opposite chiral A+(r) component is induced 
around the vortex core. At low field, the induced compo- 
nent is decayed outside of the vortex core. But at higher 
field, | A_|_ (r) | has maximum at the S-points on the bound- 
ary. The induced component A+(r) always vanishes at 
the vortex center and at the C points. At the vortex 
center (the C points), |A+(r)| recovers with the r-liner 
(r 2 -) dependence. These r-dependences are related to 
the phase winding of A+(r), which is presented in Fig. || 
(b) schematically. At the vortex center, when the dom- 
inant component A_(r) has +I(x27r) winding as shown 
in Fig. H (a) , the induced A + (r) component has opposite 
— 1 winchc^Jtis-jConsistent with the results of previous 
theories BM EIBB Since the total of the winding should 
be +1 within a unit cell, A+(r) has also +2 winding at 
the C-points, i.e., corners of the square unit cell. These 
winding structures are the same for all H range. 

The screening current |J(r)| has maximum at the scale 
of the vortex core radius, and it is decreased with ap- 
proaching the boundary of the unit cell. At the S- and C- 
points, |J(r)| = 0. By this screening current, the internal 
field is produced. It has maximum at the vortex center, 
and it is decreased outside of the core. At low field, h(r) 
has minimum at C, and monotonically increases toward 
the S-point along the boundary. But, at high field, |J(r)| 
and h(r) show anomalous behaviors at the boundary re- 
gion. The profile of h(r) has a peak at a point between 



the S and the C points along the boundary line. 

B. Vortex structure in the p+-wave domain 

Next, we consider the vortex in the p + -wave domain, 
which is shown in Fig. ^|(c) and (d). The vortex core 
shape of the dominant component shows circular shape 
at low field, as shown in the contour line of |A+(r)| in Fig. 
H (c). The induced component A_(r) becomes zero at 
the vortex center and the C-points also in this case. The 
amplitude |A_(r)| recovers with the r 2 -behavior around 
the corners C, as in the p_-wave domain case. But, the 
recovery at the vortex center shows the r 3 -behavior, in- 
stead of the r-linear. It is because A_(r) has +3 winding 
at the vortex center, as schematically presented in Fig. 
H (c). The winding —2 at the C-point is not changed. 
Compared with the p_-wave domain case, internal field 
h(r) at the vortex core is larger, since |J(r)| around the 
vortex is larger. 

At high field, the winding structure of the induced 
component A_(r) is changed around the vortex center. 
The +3 winding at the vortex center splits to a —1 wind- 
ing at the center and four +1 winding points around the 
core, as schematically shown in Fig. |Jj (d). As is seen in 
Fig. 0(d), |A_(r)| = at these +1 winding points. As 
for |J(r)| and h(r), the high field case shows the similar 
structure as in the low field case, while their strengths 
are suppressed with increasing field. 

C. In the triangular lattice configuration 

We briefly report the vortex structure in the triangular 
vortex lattice configuration. In this case, the winding 
structure at the corner of the unit cell is changed. In the 
p_-wave (p + -wave) domain case, +2 (—2) winding in the 
square lattice becomes +1 (—1) winding at the corners of 
the hexagonal unit cell in the triangular lattice, as shown 
in Fig. [| (a)-(c). In the square lattice case, the winding 
structure at the vortex center changes from +3 to — 1 on 
raising field in the p + -wave domain case [Fig. [| (c) and 
(d)], the vortex center keeps +3 winding for all H range 
in the triangular lattice configuration. 

As an example, we show the profile of the vortex struc- 
ture for the p_-wave domain at H = 0.8 in Fig. [|(d), 
which is plotted along the NN direction (VS line in Fig. 
f|(a)), boundary line (SC line) and the NNN direction 
(VC line) in the triangular vortex lattice. Compared 
with Fig. |(d) in the square lattice configuration, the 
recovery of the induced component around the C-points 
shows r-linear relation instead of the r 2 -behavior, reflect- 
ing the change of the winding structure. There appears 
anomalous field distribution at higher field in the p-- 
wave domain. But, its profile is different from that of 
the square lattice case. In the triangular lattice, h(r) has 
minimum at the S-points. However, the p+-wave domain 
case and the low field p_-wave domain case show quali- 
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FIG. 2: Vortex structure in the square vortex lattice configuration of the p_-wave domain case at low field H = 0.1 (a) and 
at high field H — 0.8 (b), and in the p+-wave domain case at H = 0.1 (c) and H = 0.4 (d). From upper panels, we plot 
the dominant pair potential |A=p(r)|, the induced component |A±(r)|, the screening current |J(r)|, and the internal field h{v), 
respectively. The left panels are the stereographic view within a unit cell region. There is a vortex at the center of the figures. 
The right panels are the profiles along the path V-S-C-V presented in Fig. pl(a). 





(d) »-l + 1 



FIG. 3: Phase winding structure of the pair potential in 
the square vortex lattice configuration. The vortex centers 
and other singularity points are, respectively, presented by 
• and o. We show the winding number for these points in 
the figures. The square area means a unit cell of the vortex 
lattice, (a) The phase of the dominant component. There is 
a winding +I(x27r) at the vortex center. Along lines V-S- 
C-V, we consider profiles of the vortex structure in Fig. ^[ 
(b) The phase of the induced pair potential A+(r) in the p~- 
wave domain case, (c) The phase of the induced pair potential 
A_(r) at low field in the p+-wave domain case, (d) The same 
as (c), but at high field. Total winding number is +3 within 
the dotted circle around the vortex core. 




FIG. 4: Phase winding structure in the triangular vortex 
lattice configuration for the dominant component (a), and 
the induced component in the p_-wave domain case (b) and 
the p+-wave domain case (c). The vortex centers and other 
singularity points are, respectively, presented by • and o. We 
show the winding number for these points in the figures. The 
hexagonal area means a unit cell of the vortex lattice, (d) 
The profile of the vortex structure in the triangular vortex 
lattice configuration of the p_-wave domain case at high field 
H = 0.8. We plot the dominant pair potential |A_(r)|, the 
induced component |A+(r)|, the screening current |J(r)| and 
the internal field h(r) along the path V-S-C-V presented in 
(a). 
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FIG. 5: Field dependence of the amplitude of the dominant 
and the induced pair potentials in the p_-wave domain case 
(a) and in thep+-wave domain case (b). We plot max|A_(r)| 
and maxjA+(r)j as a function of H. Solid (dashed) lines are 
for the square (triangular) vortex lattice configuration. 



tatively the same field distribution as that of the square 
lattice case in the profile plot. There, h(r) has minimum 
at the C points. 



D. Magnetic field dependence 

In this subsection, we investigate the continuous field 
dependence of the vortex structure. The field dependence 
of A™ ax = max[A_(r)| and A™ ax = max|A + (r)| is pre- 
sented in Fig. [|. We show it both for the square lattice 
and the triangular lattice configurations. In the p_-wave 
domain case, the induced A™ ax is large. When the am- 
plitude of the dominant A™ ax is decreased with raising 
field, the ratio of the induced component, A™ ax /A™ ax , 
increases monotonically up to H c2 . Due to this large in- 
duced component, the superconductivity in the p_-wave 
domain case can survive until higher magnetic field, giv- 
ing high H c2 . 

On the other hand, in the p+-wave domain case, the 
induced component A™ ax is small. The ratio of the in- 
duced component, A'" ax /A™ ax , decreases as a function 
of H at higher field, after increasing at low field. Since 
the ratio is reduced to zero at H — > H C 2, the p+-wave 
domain has the same H C 2 as in the isotropic s-wave pair- 
ing in the two dimensional Fermi surface case. We also 
calculate the isotropic s-wave case, which is equivalent to 
the case when we neglect the induced component of the 
chiral p-wave case. The field dependence of max|A(r)| in 
the s-wave pairing is almost the same as that of Al? ax in 
Fig. U (b) . The large amplitude of the induced compo- 
nent in the p_-wave domain is the origin of the different 
behavior of the vortex structure between the p+-wave 
domain and the p_-wave domain cases. 

The magnetic field dependence of the internal field 
distribution is shown in Fig. |^. There, we plot in- 
dependence of the maximum /i max and the minimum /i m ; n 
of h(r). The p + -wave case and the s-wave pairing case 
have similar internal field distributions. h max is the in- 
ternal field at the vortex center. Compared with the 
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FIG. 7: Zero-energy local density of states in the p_-wave 
domain case at H =0.1 (a), 0.8 (c), and in the p+-wave do- 
main case at H =0.1 (b), 0.4 (d). We plot N(E = 0,r)/N 
within a unit cell. The peak at the vortex center is truncated 
at N(E = 0,r)/iVo = 2 to show the tail structure extending 
from the vortex core. 



FIG. 6: Field dependence of the internal field distribution, 
(a) Internal field at the vortex center, /i max , in the p_-wave 
and the p+-wave domain cases. Both the square and the tri- j 
angular lattice configurations give same /i msx . (b) The min- ~s 
imum /i m in of the internal field, and the peak field h a of the 
distribution function P(h) in the p--wave domain case. Solid ^ 
(dashed) lines are for the square (triangular) vortex lattice 1 
configuration, (c) The same as (b), but in the p+-wave do- 
main case, (d) Internal field distribution function P(h) at 
H — 0.4 in the p+-wave domain case. Solid (dashed) lines are 
for the square (triangular) vortex lattice configuration, (e) " 
The same as (d), but at H = 0.8 in the p_-wave domain case. ( a ) 
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p_-wave domain case, /i max in the p + -wave domain case 
is larger at low H , and becomes smaller at higher H, be- 
cause /i max - > with approaching low H C 2- We obtain 
same h max both in the square and the triangular vortex 
lattice configurations. 

The internal field distribution function is defined as 

P(h) = {S(h-h(v))) I (15) 

with the 5-function. We can observe P{h) as the reso- 
nance line shape in the NMR or /iSR experiments. Figure 
^(d) shows P{h) at H = 0.4 in the p + -wave domain case. 
It is a typical distribution shape as in the conventional su- 
perconductor, i.e., /i m in is the internal field at the C-point 
on the unit cell boundary, and peak field h—corresponds 
to the internal field at the saddle point S.cfl The vortex 
lattice configuration affects the distance between h s and 
^min- Compared with the square lattice configuration, 

— ^min is small in the triangular lattice configuration. 
Figure ||(c) shows that this property of P(h) appears for 
all H region in the p + -wave domain case. 

Figure |(e) shows P{h) at H = 0.8 in the p_-wave 
domain case, which gives anomalous distribution of h(r) 



FIG. 8: Spectrum of the spatially averaged DOS, N(E) /N . 
(a) In the p_-wave domain case at H =0.05, 0.1, 0.4, 0.8. (b) 
In the p+-wave domain case at H =0.05, 0.1, 0.2, 0.4. The 
square and the triangular vortex lattice configurations give 
same result. The spectrum of the p+-wave domain case is 
almost the same as in the s-wave pairing case. 



at the unit cell boundary, as shown in Figs. |^(b) and 
||(d). There, the position giving h min (h s ) shifts to other 
position than the C- (S-) point on the boundary. In this 
distribution, the square and the triangular lattice con- 
figurations give similar P(h) distribution. When we see 
the field dependence of h s and /i m in in Fig. ^(b), there is 
a small difference between the triangular and the square 
lattice configurations at higher field H > 0.4. But, at 
lower field H < 0.4, there appears eminent difference in 
the distance between h s and h m i n for the two vortex lat- 
tice configurations. There P(h) has similar distribution 
as in Fig. ^|(d). It is interesting that the square vortex 
lattice configuration has lower free energy than the trian- 
gular one in the field region H > 0.4, showing anomalous 
field distribution. 
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FIG. 9: Field dependence of the zero-energy DOS N(0)/N 
in the p_-wave domain case (•), the p+-wave domain case 
(o), and the s-wave case (*). Points show numerical data, 
and lines are fitting curves by N(0) oc H a . Dashed line and 
dotted line are for a =0.74 and 0.71, respectively. Solid line 
is for N(0) oc \[H. Inset: N(Q)/No is plotted as a function of 
yH. The solid line shows the relation N(0) oc yH. 




FIG. 10: Effective pair potential for zero-energy quasiparti- 
cles. We plot |A cff (r)| = |A(0 r ,r)| at low field H = 0.1 (a) 
and at higher field H — 0.8 (b) in the p_-wave domain case, 
and at H — 0.4 in the p+-wave domain case (c). The left pan- 
els are the stereographic view within a unit cell region. There 
is a vortex at the center of the figures. The right panels are 
the profiles along the path presented in Fig. 0(a). 



V. QUASIPARTICLE STRUCTURE IN THE 
VORTEX STATE 

The LDOS is expected to be observed by the scanning 
tunneling microscopy (STM), which will experimentally 
give the detailed information of the quasiparticle struc- 
ture around vortices. We study the LDOS at zero energy, 
which dominantly contributes to the low temperature be- 
haviors. Figures [?] (a) and (b) show N(E — 0, r) at low 
field H = 0.1 in the p_-wave and the p+-wave domain 
cases, respectively. Since the gap function has full gap as 
in the s-wave pairing, the low energy states are localized 
at the vortex core. We see that the localized LDOS is 
slightly suppressed along the NN and NNN vortex direc- 
tions. It is the effect of thxjJjpi^YQrtex transfer of the 
low energy quasiparticles.r 5 H 16 H 17 li ls l'Ej 

Higher field case at H ~ 0.5H c2 are presented in Fig. 
(c) and (d). Figure (c) is for the p_-wave domain 
case at H = 0.8, and Fig. |^ (d) is for the p + -wave do- 
main case at H — 0.4. Since the inter-vortex distance 
becomes short at high field, the LDOS localized at vor- 
tex cores are overlapped each other. We see the eminent 
suppression along the NN and NNN vortex directions due 
to the inter-vortex transfer. The localized LDOS around 
the vortex core is reduced to uniform distribution when 
approaching H c2 . Then, the sharp peak at the vortex 
center survives until higher field in the p_-domain case, 
since it has higher H c2 . 

The spectrum of the spatially averaged DOS, N(E), is 
presented in Fig. ||. The full gap structure for E < Ao 
and the peak of the gap edge at E = A are gradu- 
ally smeared by low energy quasiparticles around the vor- 
tex. These behaviors of the LDOS and spectrum N(E) 
are the sarae.|ajL in the isotropic s-wave case previously 
reported.E3t§EZliia When we see the spatial structure of 
the LDOS, we do not find drastic changes by the chirality 
effect. We see qualitatively the same structure both for 
the p_-wave domain and p+-wave domain cases. 

However, when we consider the quantitative field de- 
pendence of the DOS N(0) by spatially averaging the 
LDOS at E — 0, we can see the characteristic behavior of 
the chiral p-wave pairing. Figure |^ shows the field depen- 
dence of iV(0). Numerical data are presented by points, 
and lines are fitting curves by the relation N(0) oc H a . 
The square and the triangular vortex lattice configura- 
tions give the same result. 

In the p+-wave domain case, N(0) shows similar be- 
havior as that of the s-wave pairing case. But N(0) is 
slightly smaller than that of the s-wave case at low field. 
Fitting curve is given by a = 0.74 in the s-wave case for 
K = 2.7. In the p_-wave domain case, lower field data 
are fitted by a — 0.71. But, higher field data are fitted 
by a = 0.5. To see the \/i?-behavior, we plot N{0) as a 
function of s/H in the inset. At higher field, N(0) is on 
the line iV(0) oc y/~H. This \/H- behavior at higher field is 
qualitatively fS«nsistent with the experimental data of the 
specific heat.EJ We note that this \fll is not the so-called 
Volovik effect for the vertical line node of the k 2 , — k 2 - 
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type.t§E3 In the vertical line node case, the -\Zff-behavior 
appears from low field. According to the recent direc- 
tional dependent thermal-conductivity experiment under 
parallel field,n3 there is no in-plane gap anisotropy, im- 
plying the absence of the vertical line node. 

Unless we consider the induced component of the pair 
potential, |A(0, r)| = |A_(r)||</>_(0)| gives full gap, as in 
the s-wave pairing case. Then, we can not expect the 
\/i?-bchavior of N(0). However, when we take account 
of the induced component, |A(0, r)| = |A_(r)<^>_(0) + 
A + (r)<j> + (6)\ can be small for particular direction 9, if the 
induced component |A + (r)| is large. To discuss the origin 
of the v^ff-behavior at higher field, we consider the effec- 
tive pair potential for zero-energy quasiparticles. When 
we analyze the LDOS at E = 0, the main contribution 
comes from the quasiparticle trajectory passing through 
the vortes center (i.e., line with the impact parameter 
rj_ = 0),E3 since zero-energy quasiparticles are localized 
around the vortex core. Then, the zero-energy LDOS at 
r = (x,y) dominantly consists of quasiparticles travel- 
ing along the quasiparticle trajectory with the direction 
9 r = tan _1 (y/x). They feel the effective pair potential 
A off (r) ee A(0 r ,r) - A+(r)0+(0 r ) + A_(r)0_(0 r ). 

Figure [To] shows amplitude |A c g(r)| for some typi- 
cal cases. In the p^-wave domain case, |A c g(r)| ~ 
|A_(r)| — |A + (r)|, i.e., the induced component |A + (r)| 
suppresses the effective pair potential. At low field, since 
the induced component is localized around the vortex 
core, |A c g(r)| is suppressed around the vortex core as 
shown in Fig. |o|(a). Since the amplitude of the effective 
pair potential is still large at the boundary of the unit 
cell, we expect the bound quasiparticle states around 
the vortex core. Then, the exponent a ~ 0.71 is not 
largely different from that of the s-wave case. At higher 
H, the ratio of the induced component is enhanced, and 
|A+(r)| has large amplitude at the boundary of the unit 
cell. Since the induced component |A+(r)| is large in the 
NN vortex direction [Fig. |](b)], |A e ff(r)| is largely sup- 
pressed in this direction, as shown by line VS in Fig. 
Qb). The effective pair potential is not suppressed 
at the C-points, since there is no induced component 
there. Also in the triangular vortex lattice configura- 
tion, |A c g(r)| gives the similar structure. There, |A e gf(r)| 
is eminently suppressed along six NN vortex directions. 
With increasing H, |A e gf(r)| along the NN direction is 
more suppressed, since the ratio Af ax /A max increases 
monotonically. Then, the low energy quasiparticles can 
easily transfer between NN vortices at higher field. When 
low energy quasiparticles are extended to the boundary 
of the unit cell as in the rf-wave pairing case, we expect 
the vjff-like behavior.li§E3 It is the origin of the relation 
N(0) ~ y/H at high field. 

On the other hand, in the p+-wave domain case, 
|A eff (r)| ~ |A+(r)| + |A_(r)|, i.e., the induced com- 
ponent |A_(r)| enhances the effective pair potential, as 
shown in Fig. IQ(c) . Then, low energy quasiparticles are 
bound states, as in the s-wave pairing case in all H range. 
Due to the enhancement of the effective pair potential by 



the induced component, N(0) is slightly suppressed than 
that of the s-wave case. With approaching H c2 , N(Q) is 
reduced the s-wave case's value, since the induced com- 
ponent is decreased to zero. 



VI. SUMMARY AND DISCUSSIONS 

We have investigated the field dependence of the vor- 
tex structure in chiral p-wave superconductors. We have 
shown the difference of the vortex structure for the p + - 
wave domain case and the p_-wave domain case. The 
difference comes from the structure of the induced op- 
posite chiral component. When we compare the free en- 
ergy, the p_-wave domain is the stable state, and the 
p_l_-wave domain is metastable. Then, the transition of 
the p+-wave domain to the p_-wave domain may occur. 
We expect different vortex lattice configuration for the 
domains, i.e., square- like lattice in the p_-wave domain 
and triangular lattice in the p+-wave domain. 

The phase winding structure of the induced component 
of the pair potential is different depending on the chiral- 
ity. In the p+-wave domain, the amplitude of the induced 
component is small and reduced to zero near H C 2 . Then, 
the vortex structure is similar to that of the isotropic 
s-wave pairing case. And H c i is same as in the s-wave 
pairing case in the two-dimensional Fermi surface. In the 
p__-wave domain case, the opposite chiral component is 
largely induced. Then, the superconductivity can survive 
until high field, giving high H c %. The induced component 
produces the characteristic vortex structure in the chiral 
p-wave superconductors, such as an anomalous internal 
field distribution. 

The LDOS structure shows that low energy quasi- 
particles are bound states around the vortex core, and 
there are some inter- vortex transfers. At higher field, 
the bound states are overlapped between neighbor vor- 
tices. When we quantitatively consider the field depen- 
dence of the zero energy DOS iV(0) oc H a , we obtain the 
effect of the chiral p-wave superconductivity. The sta- 
ble p_-wave domain case shows -\/ff-behavior at higher 
field. It is because the effective pair potential for zero 
energy quasiparticles are suppressed along the NN vor- 
tex direction by the induced opposite chiral component 
of the pair potential. At low field, the suppression by the 
induced component is restricted in the vortex core region, 
the low energy quasiparticles are still bound around the 
vortex core. Then the exponent a is near the value for 
the s-wave pairing at low field. 

The superconducting state in Sr2Ru04 is suggested 
to be the chiral p-wave pairing. If we can experimen- 
tally observe the domain structure of the p+-wave and 
the p_-wave pairing regions, it becomes firm evidence 
for the chiral p-wave superconductivity. In this observa- 
tion, the information of the vortex structure difference 
for the p±-wave domains is helpful to analyze the chiral- 
ity p± of each domain. For example, the internal field 
distribution and the stable vortex lattice configuration 
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may be different depending on the chirality of the do- 
main. The specific heat measurement on Sr 2 Ru04 re- 
ports that N(0) (x \IS at high field, while it deviates 
from \J~~H at low field.t£l It is qualitatively consistent with 
our results. However, when we analyze the experimental 
data on Sr2Ru04, there are some factors to quantita- 
tively modify our results on a simple isotropic system, 
such as the possibility of the line node along the basal 
plane direction, the orbital dependence and anisotmpv nf 
the Fermi surface and the gap functionsEJSE3EJoE3'Lll 
The study on these additional effects remains in future 
problems. 
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(f>(—9) = e la (f>*(9) for the pairing function of the domi- 
nant component of the pair potential. Then, considering 
the transformation Sr and 9 — > — 9 in Eqs. (0)-(§), we 
obtain the following relations of the quasiclassical Green's 
functions, 

/K,-0,Sr) = -e itt r( W „,0,r), 
/t(w n ,-0,Sr) = -e- iQ /t*K,0,r), 
g(u n ,-0,Si)=g*(u> n ,6,T)- (A2) 

Next, we consider the V'-rotation around the vortex 
center at r = 0, i.e., R^,r — (a;cosV> — y sin ip, x sin?/) + 
ycosip). Generally, vortex lattice has the symmetry for 
the rotation ip = tt. And further, the square (triangular) 
vortex lattice has the symmetry for the rotation ip = tt/2 
(ifj = 7r/3). Under these rotations, A(9 + ip,R^,r) = 
e i(a -V>)A(0, r) in the symmetric gauge. The factor e' a 
comes from <f>(6 + if) = e la (f)(0). Then, considering the 
translation R^r and 6 — > 9 + ip in Eqs. (0)-(||), we obtain 



APPENDIX A: SYMMETRY RELATION 

When one of the vortex center locates at r = 0, there 
is a relation A(9,r) = —A(0, — r). Then, considering the 
transformation r — > — r in the Eilenberger equations ([!])- 
(|J) , we obtain the following relations of the quasiclassical 
Green's functions, 

/K,0,-r) = -/t*K,0,r), 
/tK,9,-r) = -/K,«,r), 
g(u n ,0,-r)=g*(u* n ,0,r). (Al) 

Then, in the calculation of the Matsubara frequency Lo n 
or E = 0, it is enough to solve the Eilenberger equations 
in half area of a unit cell. 

When the vortex lattice is symmetric under the reflec- 
tion at the x axis, i.e., Sr = (x, —y), there is a relation 
A(-9,Sr) = -e ia A*(9,r). The factor e iQ comes from 



/K, 9 + 4>, iv) = ^"'"/K, 0, r), 

P(oj n , 9 + ^ R 4 ,r) = e-^ a '-^p(iu n , 9, r), 
g(u n , 9 + ip, R^r) = g(u n , 9, r). (A3) 

In the pair potential A(0,r) = A_(r)0_(0) + 
A + (r)4> + (9), the pairing function <f>±(9) of the induced 
component may produce different phase factor from that 
of the dominant component tf>^ (9) in the transformation 
(f>(9 + ip) — e IQ (f>(9). Then, the phase of the induced 
component A±(r) should be changed so as to cancel the 
difference of the phase factors in the rotational transfor- 
mation. This is the origin of the different phase winding 
of the induced component in Figs. [5] and ^. 

Using these symmetry relations, it is enough to solve 
the Eilenberger equations (0)-(||) for < 9 < n/4 
(0 < 9 < 7r/6) in the square (triangular) vortex lattice 
configuration. 
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